Bulk and surface magnetoinductive breathers in binary metamaterials 
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We study theoretically the existence of bulk and surface discrete breathers in a one-dimensional 
magnetic metamaterial comprised of a periodic binary array of split-ring resonators. The two types 
of resonators differ in the size of their slits and this leads to different resonant frequencies. In the 
framework of the rotating-wave approximation (RWA) we construct several types of breather exci- 
tations for both the energy-conserved and the dissipative-driven systems by continuation of trivial 
breather solutions from the anticontinuous limit to finite couplings. Numerically-exact computa- 
tions that integrate the full model equations confirm the quality of the RWA results. Moreover, it 
is demonstrated that discrete breathers can spontaneously appear in the dissipative-driven system 
as a results of a fundamental instability. 
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I. INTRODUCTION 



Discrete breathers (DBs) or intrinsic localized modes, 
are time-periodic and spatially localized excitations that 
may be produced generically in discrete arrays of weakly 
coupled nonlinear elements [H, [1] ■ A large body of the- 
oretical work has produced means of precise numerical 
analysis of coupled oscillator systems with breathers both 
in the Hamiltonian as well as dissipative regimes 0, 
S H) Breathers have been experimentally observed 
in several diverse systems, including optical waveguides 
systems [Ij, solid-state systems jj antiferromag- 
netic chains [ll| , Josephson junction arrays [l^ , and mi- 
cromechanical oscillators [ij, among others. Dis- 
crete breathers can be generated spontaneously in a lat- 
tice either through stochastic mechanisms [lal or by 
purely deterministic mechanisms [TtI [TsI . Il9l . l2Cll | in a 
process by which energy, initially evenly distributed in 
a nonlinear lattice, can be localized into large ampli- 
tude nonlinear excitations. Indeed, it has been demon- 
strated experimentally that the standard modula- 
tional instability (MI) mechanism in dissipative systems 
driven by an alternating term can initiate that process 
by the formation of low-amplitude breathers. The en- 
ergy exchange between those low-amplitude DBs favors 
the higher-amplitude ones, resulting eventually in the for- 
mation of a few high-amplitude DBs. 

A few years ago a whole new class of artificially struc- 
tured materials, referred to as metamaterials, was dis- 
covered; the latter are comprised of discrete elements 
and exhibit electromagnetic properties not available in 
naturally occurring materials. A subclass of those meta- 
materials, the magnetic metamaterials (MMs), exhibit 
significant magnetic properties and sometimes even neg- 
ative magnetic resp onse up to Terahertz (THz) and opti- 
cal frequencies [2ll [22| . The most common realization of 
a MM is composed of periodically-arranged electrically 



small sub-wavelength particles, referred to as split-ring 
resonators (SRRs) [H, [2^. In its simplest form, each of 
those resonators is just a highly conducting metallic ring 
with one slit. The MM thus built can become nonlinear 
either by the insertion of a nonlinear dielectric [25| or a 
nonlinear electronic component (e.g., a varactor diode) 
[26I \2T\ in the slit of each SRR, resulting in voltage- 
dependent SRR capacitance. In microwave frequencies, 
such a MM has been realized [2^ and is dynamically 
tunable by varying the input power. The combination 
of nonlinearity and the discreteness that is inherent in 
those metamaterials, makes possible the generation of 
nonlinear excitations in the form of DBs. The existence 
and stability of DBs in nonlinear SRR-based MM mod- 
els, that are localized either in the bulk [2^, |3^ or at 
the surface [3l|, [s^l of the MM, have been demonstrated 
numerically. Moreover, domain walls [ssj and envelope 
solitons [34] may also be excited in such systems. The 
surface DBs in MMs are very similar with the surface 
modes observed in discrete waveguide arrays [H, [s^ . 

Recently, a novel MM comprised of two types of SRRs 
was investigated theoretically [37|, and it was demon- 
strated that in the nonlinear regime such binary MMs 
are perfectly suited for the observation of phase-matched 
parametric interaction and enhanced second harmonic 
generation (SHG). In the present work we extend the 
previous studies on DB generation in model SRR-based 
MMs to the case of a one-dimensional (ID) binary MM 
with on-site cubic nonlinearity. In practice, a binary MM 
can be constructed in many different ways, by changing 
for example one or more of the material and/or the ge- 
ometrical parameters of the SRRs belonging to one type 
(say type a), with respect to the same parameters of the 
SRRs belonging to the other type (say type b). Here we 
allow for different slit-widths for the two types of SRRs, 
which makes their resonance frequencies differ by a fac- 
tor that we call resonance frequency mismatch (RFM). 
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FIG. 1: One-dimensional binary array of split-ring resonators. 



SRR radius, h is the diameter of the metal wire, d is the 
slit width of the SRR, and p the (material-dependent) 
SRR resistivity. Neighboring SRRs in an array are mag- 
netically coupled through their mutual inductance M, 
that is approximately given by 



The considered binary MM is formed by type a SRRs 
at the even-numbered sites and type b SRRs at the odd- 
numbered sites of a periodic ID array. In the next section 
we give the model equations that describe the dynamics 
of the binary MM and we obtain the corresponding linear 
dispersion relation for magnetoinductive waves in that 
medium [s^. In Sections III and IV we construct, 
using the rotating wave approximation (RWA), several 
types of Hamiltonian and dissipative breathers (DDBs), 
respectively. The dynamic stability of those DBs is dis- 
cussed in Section V where the full model equations are 
integrated numerically. In most of the investigated cases 
the numerics confirm the quality of the RWA results. 
Moreover, in that Section we also demonstrate the spon- 
taneous generation of DDBs induced by MI, using fre- 
quency chirping of the driving field. That procedure has 
been used in actual experiments for the generation of 
high-amplitude DDBs in micromechanical cantilever os- 
cillator arrays [13|, and perhaps it could be also used in 
similar experiments involving nonlinear binary MMs. We 
finish in section VI with the conclusions. 



II. MODEL BINARY METAMATERIAL AND 
LINEAR DISPERSION 

Consider a ID SRR-based MM comprised of nonlin- 
ear units shown schematically in Fig. 1. Each non- 
linear SRR in the array can be mapped to a non- 
linear resistor-inductor-capacitor (RLC) circuit featur- 
ing self- inductance L, ohmic resistance i?, and nonlin- 
ear (voltage-dependent) capacitance C(|Ep) cx e(|Egp), 
where e is the field-dependent permittivity of the infill- 
ing dielectric, E is the electric field, and Eg is the electric 
field induced along the SRR slit. We assume that the lat- 
ter originates from an alternating magnetic field that is 
applied to the MM perpendicularly to the SRR planes, 
and it is proportional to the voltage U across the slit. 

Let us for the moment ignore the nonlinearity and set 
C = Ci, with Ci is the linear capacitance that is built up 
across the slit. Just like an RLC resonator, the SRRs 
exhibit an inductive-capacitive resonance at frequency 
ujR ~ Ci (for i? ~ 0, implying low Ohmic losses). 

For a circular SRR with rectangular cross-section, the 
parameters L, R, and C/ of the equivalent RLC circuit 
can be estimated from the relations [21] 
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where Eq ^-nd fio ^-re the permittivity and permeability 
in vacuum, respectively, e; is the linear relative dielectric 
permittivity of the infilling dielectric, r is the average 
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where D is their center-to-center distance, and , rj, are 
their average radii. In order to construct a binary array, 
we have to change one or more of the material and/or 
geometrical parameters of the SRRs that are going to be 
of one type, with respect to the same parameters of the 
SRRs that are going to be of the other type. As it can 
be observed from Eqs. ^ and ^ 

• A change in the SRR radius r affects L, M and R. 

• A change in h affects L, Ci and R. 

• A change in e; affects Ci which in turn, implies a change 
in the nonlinear response. 

• A change in d affects C/ and slightly R and L. 

• A change in resistivity p affects R. 

Obviously there many possibilities for constructing two 
types of SRRs and consequently a binary array. Here 
we make the relatively simple choice to create two types 
of SRRs by considering different slit-widths, i.e., da for 
type a and db for type b. Thus, the linear capacitances 
of type a and b SRRs become respectively Cq and Cb, 
resulting in different resonance frequencies uja — ^/VLCa 
and ojft = 1 / \JLCb- 

Now let us return to the nonlinear problem, and as- 
sume that the slits of all the SRRs in the array are filled 
with a Kerr-type dielectric. Then, the charge Q,,. accu- 
mulated in the capacitor of the nth SRR is [40| 



(3) 



where C„ — Ca (C„ = Cb) for SRRs at even- (odd- 
) numbered sites of the array, and x — a/(3e;) is the 
dimensionless nonlinearity coefficient, with a = +1 {a = 
— 1) for a self-focusing (self-defocusing) dielectric. The 
above equation leads to the following approximate form 
for the voltage Un across the slit of the nth SRR 
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where is a characteristic (large) voltage. Then, the 
coupled equations describing the charge dynamics in the 
nonlinear MM that is placed in an alternating magnetic 
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field are generally written as 
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where k is the normalized wavenumber, into the lin- 
earized Eqs. © and (jlOp . without the dissipative and 
the external driving terms (7 = and eg = 0) 
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where we assume that i?,„ — R, Lm — L, and 

F{t)=£osm{Lut). (7) 

The above equation gives the electromotive force (emf) 
of amplitude £q and frequency uj that is excited in each 
SRR due to the action of the field. Let us define the 
quantities 



(8) 

With the above definitions, Eqs. ([5]) and © can be writ- 
ten in normalized form as 

^[A(72n-l + q2n + M2n+l\ + <^'?2n - X^^ ^In = 

-7^ + £osin(f7r), (9) 

j2 3 
" r\ I I \ 1 I '?2n+l 92)1+1 
^[Aq2n+1 +qn + M2n+2\ H ^ 

-7— -KeoSm(S2T), (10 

where A = M/L is the diniensionless coupling param- 
eter, Q, — oj/ljq is the dimcnsionless driving frequency, 
5 = Wa/uji, is the RFM ratio, eo — £q/Uc^ and n is an 
integer. From Eqs. (O and pH)) we can see that the 
change in the linear capacitances also affects the nonlin- 
ear terms, and that, actually, the latter are affected much 
more than the linear ones. The changes that are caused 
to the terms proportional to R and L are of higher order 
and thus they are neglected. As the RFM 5 increases, 
the resonance frequency as well as the nonlinear term of 
the even-numbered site SRRs increase, while at the same 
time the resonance frequency and the nonlinear term of 
the odd-numbered site SRRs decrease. The inductive 
coupling parameter A can be either positive or negative 
depending on whether the array geometry is of the "ax- 
ial" or "planar" type [soj . 

Such MMs and other systems with magnetically cou- 
pled elements, support in the low-amplitude (linear) limit 
a new kind of waves, the magnetoinductive waves (38l.l39|. 
[3l,[3§|. In the present case, the frequency dispersion for 
linear magnetoinductive waves is obtained by substitut- 
ing 
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The dispersion curves for a particular choice of RFM 
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FIG. 2: 

and A = 



(Color online) Linear dispersion relation for S = 0.8 
-0.2. 



S and coupling parameter A are shown in Fig. 2. Even 
though we can think of each linear SRR as an RLC cir- 
cuit, i.e., an electromagnetic oscillator, we notice that 
the dispersion curves for the coupled array do not con- 
tain any acoustic-like branch; rather both curves are of 
the "optical" type: lim^^o w{k) ^ 0. This is due to the 
particular (inductive) nature of the coupling between the 
SRRs. There are now three regions where we can look 
for breathers in the nonlinear case: Above and below the 
bands, and in the Bragg gap (BG) in between. In the 
absence of RFM {6 = 1) we recover the single band for 
the 'monoatomic' MM [301 • As S diverges from unity the 
BG increases, pushing the f2-f branch upwards and the 
51 _ branch downwards. 

Eqs. ^ and PH)) can be written conveniently in the 
compact form 
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[Aq„_i + qn + Xqn+i] + t^„gn - XW„g; 
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dqn 
dr 
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where = (5 {ujf^ = 1/6) for even (odd) n. Without 
dissipation and external driving, the earlier equation can 
be obtained from the Hamiltonian Ti. = 7i„ , where 
the discrete Hamiltonian density 7i„ is given by 
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The last term on the right-hand side of the earlier equa- 
tion is the nonlinear on-site potential which is given by 



Vn = V{qn) = i(w„g„)^ 



1 - \x^l{^nq,if 
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The Hamiltonian 7i is actually the conserved energy of 
the lossless system in the absence of any driving terms 
system. For the dissipative system, H. is also useful since 
its time-average per period gives correctly the average 
energy per period for that system. 



III. HAMILTONIAN BREATHERS IN THE 
ROTATING- WAVE APPROXIMATION 

A standard method of DB construction in Hamilto- 
nian systems, that gives numerically exact results up 
to arbitrary precision, uses the Newton's method 0, 
which has been applied successfully for DB generation 
in MMs [H, [lOl- In this Section, we use the RWA 
method that keeps a simple physical picture and more- 
over can produce quite accurate results. According to the 
simplest version of that method, one looks for station- 
ary solutions of the system that are separable with an 
assumed time-dependence (e.g., sinusoidal) of the form 
<ln{T) = qn sin(r2r). Direct substitution of QniT) into Eq. 
(fTS]) with the approximation sin(a;)^ « (3/4)sin(a;) gives 
an algebraic system of nonlinear equations for the q„s 
that reads 



rt'^{Xqn+i + q-n + A(7„-i) + uj^q 



In the anticontinuous limit (A 
has the solutions 



, - (3/4)xc.^g3 = (|16) 
0) the earlier equation 



Qn = 



(3/4)x^r 
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According to Eqs. ^T7\ . we have the following interesting 
possible scenarios: 

(i) a > 0. Then, ql > for aU n, if < Mm{S, 1/6}. 

(ii) a < 0. Then q^^ > for all n provided Jl^ > 
Max{(5,l/(5}. 

(iii) In the intermediate case where Min{5,l/S} < 
fi^ < Max{(5, 1/(5}, what happens is that q2n > and 
q2n+i = 0, or the converse, depending upon the sign of 
a. 

The RWA method can be used for the construction of 
DBs both on the 'surface' and the bulk of the energy- 
conserved binary MM. For a ID MM, a surface localized 
DB obviously corresponds to an edge state, i.e., a state 
with maximum amplitude at either of the two ends of the 
array. A bulk DB, on the other hand, is meant to be a DB 
whose maximum amplitude is far from the end-points of 
the array. However, the procedure of obtaining DBs by 
the RWA method, both at the surface or in the bulk, pro- 
ceeds in the same way. The first step is to set up a trivial 
DB. We first choose its central site, i.e., the site where 
the DB shall have its maximum amplitude. Suppose that 
the central site is taken at n — ub where the 'coordinate' 
qn = quB 8.nd set all the g„ for n ^ ub equal to zero. 
The value of qng is calculated from Eq. P7)l . with an 
appropriately chosen 17. That solution is subsequently 
continued for finite couplings up to a maximum value 
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FIG. 3: (Color online) Typical Hamiltonian surface breather 
profiles in a magnetoinductive binary chain for S — 0.9, A = 
0.1, Q = 0.77, X = 1/6, obtained by the RWA. 



A — Xmax where DBs cease to exist. Usually we consider 
a DB to be localized around the site where it exhibits 
its maximum amplitude. For a finite array, the bound- 
ary conditions that should be imposed to the dynamic 
equations resulting from the Hamiltonian (|14p should be 
specified. For DBs excited in the bulk one may use either 
periodic or open-ended boundary conditions, since DBs 
are highly localized entities and are not affected by the 
boundaries. However, for surface DBs the termination of 
the structure should be taken into account, and for that 
purpose we should use open-ended boundary conditions. 
In the following, we always use that type of boundary 
conditions, i.e., qo = 0, qN+i ~ 0, where N is the total 
number of SRRs in the binary array. 

Surface breathers.- Typical surface- localized, single- 
site DBs, that are generally very similar to the ones 
examined for a discrete nonlinear Schrodinger (DNLS) 
model for a semi-infinite binary waveguide array |4l| , are 
displayed in Fig. 3. The unstaggered modes shown there 
originate in the lower gap region < < Mzn[(5, 1/(5]. 
There are also staggered modes (not shown) originating 
from the Bragg gap re gion that constitute magnetoinduc- 
tive Tamm states [3ll |32|| . From the surface DBs shown 
in Fig. 3, only one of them (Fig. 3a) corresponds to a 
truly surface state, since it is localized exactly at the left 
end of the array (n = 1). The next two (Figs. 3b and 
3c) can also be characterized as surface DBs, since they 
are localized very close to the surface (n = 2 and n = 3, 
respectively), but actually they are cross-over states be- 
tween surface and bulk DBs. Since the DBs shown here 
are highly localized, they obtain their bulk form within a 
distance of only a few sites from the surface, so that the 
DB shown in Fig. 3d (localized at rt = 4) can be consid- 
ered as a bulk DB. With the appropriate choice of the 
initial conditions we may also construct multi-site surface 
DBs such as those shown in Fig. 4, that remind four-site 
antisymmetric DB excitations. Of course, close to the 
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FIG. 4: (Color online) Typical Hamiltonian surface multi- 
site breather profiles in a magnetoinductive binary chain for 
S = 0.9, A = -0.1, n = 0.77, X = 1/6, obtained by the RWA. 
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FIG. 5: (Color online) Typical Hamiltonian bulk symmetric 
breather profiles in a magnetoinductive binary chain for 5 = 
0.8, X = -1/6, w = 1.5, and (a) A = +0.1; (b) A = -0.1. 
These profiles are also obtained by the RWA. 



surface, we cannot obtain exact antisymmetric DBs. 

Bulk breathers.- Typical bulk Hamiltonian DB profiles 
obtained with the RWA method are shown in Figs. 5 
and 6. In Fig. 5, the two single-site symmetric DBs dif- 
fer in that the first one (Fig. 5(a)) is staggered, while 
the other one (Fig. 5(b)) is unstaggered. The stag- 
gered/unstaggered character of those Hamiltonian DBs 
depends on the sign of the product of the coupling pa- 
rameter and the nonlinearity parameter, a = sgn{aX). 
For fT > 0, that implies either a = -1-1 and A > or 
a = —1 and A < 0, the excited DBs are unstaggered. On 
the other hand, for ct < 0, that implies either a — +1 
and A<OorQ; = — 1 and A > 0, the excited DBs are 
staggered. In that figure the DB frequency ^Ib = '^'^/Tb, 
with Tb the DB period, was chosen to ensure that the 
DB amplitude at all sites is nonzero in the anticontinuous 
limit. In Fig. 6, that shows two-site antisymmetric bulk 
DBs, the nonlinearity parameter is negative {a = —1) 
implying that in the anticontinuous limit, only the even 
sites of the array can have a nonzero amplitude. Again, 
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FIG. 6: (Color online) Typical Hamiltonian bulk antisym- 
metric breather profiles in a magnetoinductive binary chain 
for S = 0.5, X = -1/6, Lo = 0.8, and (a) A = -0.2; (b) 
A — +0.2. These profiles are also obtained by the RWA. 



with appropriate choice of the initial conditions and by 
changing the sign of the coupling parameter A and/or the 
nonlinearity parameter a we may also construct multi- 
site bulk DBs with different symmetry. It should be also 
possible to generate a large variety of surface and bulk 
Hamiltonian DBs in two dimensional binary arrays, just 
like in planar arrays of identical SRRs |30j . Increased di- 
mensionality offers more possibilities for generating dif- 
ferent DB types. 



IV. DISSIPATIVE BREATHERS IN THE 
ROTATING WAVE APPROXIMATION 



We consider the more realistic case of dissipative dis- 
crete breathers (DDBs), that can be excited either on 
the 'surface' or the bulk of a ID binary MM. In typi- 
cal experiments that involve MMs, the metamaterial is 
driven by an applied electromagnetic field of appropriate 
polarization. In ID there are two possible geometries for 
the arrangement of the SRRs (3Q] ; the planar geometry 
(as shown in Fig. 1), for which the coupling parameter 
is negative, and the axial geometry for which the cou- 
pling parameter is positive. That polarization can be 
chosen such that, for example, the magnetic component 
of the field is perpendicular to the SRR planes, while the 
electric component is parallel to the sides of the SRRs 
that do not have a slit. This choice simplifies physically 
the situation, since only the magnetic field excites an 
emf in the SRRs. Thus, in the equivalent circuit pic- 
ture, a binary SRR-based MM in an electromagnetic field 
can be described by an array of nonlinear RLC circuits 
driven by an alternating emf that are coupled through 
their mutual inductances. The losses of the SRRs can 
be described in this picture by an equivalent resistance. 
That effective resistance R may actually describe both 
Ohmic losses of the SRR as well as radiative losses, if 
these are relatively low [sj]. Under those assumptions, 
the dynamics of the charges g„(T), n — 1,2,. ...,N, is 
given by Eq. ^3]) . In the framework of the RWA method, 
we look for stationary solutions of that equation in the 
form (7„(t) = g„ sin(rir + 0„). By direct substitution of 
Qn (t) into Eq. (fT51) and by making the RWA replacement 
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FIG. 7: The intersection of P{qn) (solid line) and Eq (dashed 
hne) provides the nonlinear attractor(s) for a single SRR. Case 
displayed here corresponds to 5 — l,x = 1/6,7 — 0.02, = 
0.92, eo ~ 0.04. The black (red) circle represents a stable 
(unstable) attractor. 



sin{QT + (pn)^ « (3/4) sin(17r + 0„), where g„ thereafter 
denotes the time-independent DB amplitude at the nth 
site and its phase, we find that the DB amplitudes 
and phases at each site n satisfy the relations 



+1 



by varying a parameter in that four-dimensional param- 
eter space, two of these solutions may disappear through 
a pitchfork bifurcation, leaving behind only one single at- 
tractor. The boundary in parameter space that separates 
those two cases can be found implicitly by computing the 
values of qn denoted by q^, for which dP(qn)/dqn = 
(i.e., the values of qn that correspond to the local extrema 
of P(g„)). Obviously, if qn — qn corresponds to a local 
minimum (maximum) then for Piq^) > {P{ln) < ^o) 
there is only one attractor left. Thus, the earlier inequal- 
ities determine different regions in the parameter space 
where, depending on the values of the parameters, we 
may have either one or three attractors. Thus, we may 
distinguish on such a diagram two different 'phases' that 
correspond to either three or one real solutions. A typ- 
ical example is shown in the reduced eg — f2 parameter 
space in Fig. 8, for ^ = 1 (no RFM) and opposite val- 
ues of X, while the values of the driving field strength 
eo and the driving frequency f2 are varying. There we 
see clearly the areas where there are either three (inside 
colored area) or one (outside colored area) attractor (s). 
Another example, where the RFM changes from 5 = 0.5 
to 1.5 is shown in Fig. 9 for 7 = 0.1 and positive x- 
There we observe that the colored area, corresponding 
to three attractors, expands with increasing 5. The 
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where n = 1,2, . . . , N and qo = qN+i = 0. The inclusion 
of dissipation and external driving alters significatively 
the possible DB modes that the binary SRR system can 
support. The dissipative DBs possess the character of 
an attractor for initial conditions in the corresponding 
basin of attraction, and they may appear as a result 
of power balance between the incoming power and the 
intrinsic power loss [3, Q Dissipative DB excitations in 
SRR-based MMs are of great importance since they alter 
locally the magnetic response of the system from diamag- 
netic to paramagnetic or vice versa [l^, [13, [13] . 

In the anticontinuous limit Eqs. (fTS)) and become 
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FIG. 8: (Color online) 'Phase diagram' in the reduced param- 
eter space £0 — for a single driven-damped SRR oscillator 
showing the regions with different number of attractors, for 
5 = 1, 7 = 0.1, and x = +1/6 (left); x = -1/6 (right). Inside 
(outside) the colored region we have three (one) attractors. 
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where we kept the subscript n to distinguish between os- 
cillators located either at an odd-numbered site {n =odd 
integer) or even- numbered site {n =even integer). The 
polynomial P{qn) is cubic in q^ for general values of the 
parameters i5, x, and 7, with P(0) = 0. Thus, there can 
be at most three real roots that correspond to attractors 
of the SRR oscillator (see Fig. 7), from which two are 
stable (unstable) and one is unstable (stable). However, 



FIG. 9: (Color online) Same as in Fig. 8 for a single driven- 
damped SRR oscillator with 7 = 0.1, x — +1/6, and (from 
left to right) 5 = 0.5; 5 = 0.75; 5 = 1; 5 = 1.5. 

values for the stable attractors predicted by the RWA 
are in excellent agreement with those obtained through 
dynamical evolution of the charge in a single SRR os- 
cillator. For instance, for the parameter set 5 = 0.8, 
n = 0.92, 7 0.01, eo = 0.04, x = +1/6, and for an 
even- numbered site n, Eq. (|20p predicts a single attractor 
at qf = ±0.582163, while for an odd-numbered site it pre- 
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diets three attractors at q° = ±1.23531, ^2 = ±1-33055 
and ^3 = ±0.0996811, of which the latter two are stable. 
The unstable attractor at Qi is not reachable through 
simple numerical integration of the dynamical equation. 
For the stable attractors, we have checked with direct nu- 
merical integration that their values are practically the 
same with those obtained from the RWA approach. In 
general, the presence of dissipation and driving severely 
limit the possible spatial profile of the breathers. The 
structures tend now to be either stron gly localized ones, 
or rather extended, like domain walls j27l . [29j . The sit- 
uation is similar for DDBs in the bulk (Figs. 10 and 11). 
As soon as 6 deviates from unity, that is, when we are 
dealing with a bona fide binary chain, the area in phase 
space with two stable attractors reduce (increase) as 5 
decreases (increases) from unity, if the site chosen is an 
even-numbered one, as can be seen in Fig. 8. The oppo- 
site behavior occurs for an odd-numbered site. 



In order to illustrate how the RWA method works in 
this case, we calculate some of typical surface DDBs. The 
calculation of bulk DDBs proceeds in the same way, by 
simply choosing the central site of the corresponding triv- 
ial DB that is located at n = ub somewhere in the bulk. 
For a given value of the RFM, we first determine first 
the attractors available for each single SRR oscillator, 
which is located either at odd- or even-numbered site. 
Then, we set up a trivial surface DB which is subse- 
quently continued for finite values of A. The continua- 
tion procedure proceeds in exactly the same way as that 
for Hamiltonian DBs, except that the relevant equations 
are now Eqs. (fTS]) and (fT9|) . Thus, we can obtain sev- 
eral types of surface DDBs for an interval of A up to a 
maximum, i.e., up to A = )^max- For example, for the 
parameter set x = +1/6,7 = 0.01, il = 0.5, 5 = 2, and 
So = 0.04, we have stable attractors Qi = ±4.06719 and 
^2 = ±0.160225 at odd-numbered sites, and qf = ±1.334 
and q = ±0.0228639 at even-numbered sites. We can 
set up a trivial surface DB localized at hb = 1 as 
qi = 4.06719, q2n = 0.0228639 and g2«-i = 0.160225 
(n > 1). For a trivial surface DB localized at ub = 2 
we may choose q2 = 1.334. q2n-i — —0.160225 and 
q2n = —0.0228639 {n > 1). Or, for a trivial surface 
DB localized at ub = 3 we may choose q2n = —0.160225, 
92n-i = -0.160225 {n ^ 2), and = 4.06719. Contin- 
uation of those trivial DDBs up to A = 0.025 gives the 
surface DDB profiles shown in Fig. 10. Of course there 
also other trivial DDB profiles that we could choose. Sim- 
ilar bulk DDBs can be obtained from the trivial DDBs 
given above only by changing ub to a value relatively far 
from the end-points. An illustrative example of a bulk 
DDB localized at an odd-numbered site is shown in the 
left panel of Fig. 11, while in the middle and right panels 
of Fig. 11 are shown two multi-site DDBs also localized 
at odd-numbered sites. The latter two DDBs have been 
obtained with appropriate choice of a trivial DDB profile. 
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FIG. 10; (Color online) Dissipative surface breather profiles 
for (5 = 2, n = 0.5, A = 0.025, x = +l/6i, 7 = 0.01, eo = 0.04, 
that are localized at n = 1 (left panel); n = 2 (middle panel); 
n = 3 (right panel). 
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FIG. 11: (Color online) Dissipative bulk breather profiles for 
S = 2, n = 0.5, A = 0.025, x = +1/6, 7 = 0.01, eo = 0.04, 
that are locafized at n = 1 (left panel); n = 2 (middle panel); 
n = 3 (right panel) . 



V. NUMERICALLY EXACT CALCULATIONS 

In this Section we construct several types of both 
energy-conserved and dissipative DBs which are local- 
ized either in the bulk or at the surface, using standard 
numerical algorithms 0, Moreover, in the case of a 
dissipative binary MM we also generate DDB excitations 
through a procedure that can be used for parameter val- 
ues where the homogeneous solution is modulationally 
unstable [l3|- In other words, we exploit MI to initiate 
spontaneous localization of energy in the binary array 

Hamiltonian breathers.- For the Hamiltonian binary 
MM, DBs can be constructed from the anticontinuous 
limit of Eqs. with eo = and 7 = where all the 
SRRs are decoupled [l^, [13] • Using Newton's method 
we have constructed several types of Hamiltonian, nu- 
merically exact DBs for the ID binary MM, for different 
parameter sets. The obtained Hamiltonian DB profiles 
are in excellent agreement with those obtained with the 
RWA method. 

Dissipative Breathers.- In order to generate DDBs we 
start from the anticontinuous limit of Eq. (|13p . where 
dissipation and driving are included. We identify stable 
attractors of each SRR oscillator, that is either located 
at odd- or even- numbered sites. For constructing a triv- 
ial DDB profile we need to find, for at least one of the 
two types of oscillators, two different amplitude stable 
attractors. For example, for the parameter set 6 = 2.0, 
n = 0.5, 7 = 0.01, £o = 0.04, X = +1/6, we obtain stable 
attractors q1 = ±4.06719 and q| = ±0.160225 at odd- 
numbered sites, and qf = ±1.334 and q" = ±0.0228639 
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FIG. 12: Dissipative breather profiles at maximum amplitude 
for several values of the coupling parameter as shown on the 
figure which are localized at n = 1 and n = 3. The parameters 
are O = 0.5, x = +1/6, 7 = 0.01, eo = 0.04, and 5 = 2. 



at even- numbered sites. Those values are practically 
the same with those obtained with the RWA method. 
We set up a trivial surface DB localized at = 1 as 
qi = 4.06719, (j2„ = 0.0228639 and g2«-i = 0.160225 
[n > 1). For a trivial surface DB localized at = 3 
we may choose qs = 4.06719, q2n = 0.0228639, and 
<Z2n-i — —0.160225 {n ^ 2). Continuation of those triv- 
ial DDBs gives surface DDB profiles up to \max — 0.19. 
Typical profiles for several values of A, both for DDBs 
localized at n = 1 and n = 3, are shown in Fig. 12. 

Another example is given for the parameter set 6 — 0.8, 
n = 0.92, 7 = 0.01, £0 = 0.04, x = +1/6, where the RWA 
method predicts a single attractor at qf = ±0.582163 at 
even-numbered site oscillators, while it predicts three at- 
tractors at qf = ±1.23531, q^ = ±1.33055 and q^ = 
±0.0996811, of which the latter two are stable, for an 
odd-numbered site oscillators. These values are also in 
agreement with those obtained by direct integration of 
the single SRR oscillators. We set up a trivial surface 
DDB localized at ns = 1 as gi = 1.33055, q2n = 0.582163 
and q2n-i — 0.0996811 (n > 1), and continue it up 
to Xjnax ^ 0.07 where DDBs cease to exist. Typical 
DDB profiles are shown in Fig. 13 for several values of 
A shown on the figure. A profile for A — 0.072 which 
is greater that Xmax, where the homogeneous solution is 
restored, is also shown in Fig. 13d. The frequency of the 
DDBs shown here is the same with that of the driver, 
i.e., Qb = 2it/Tb — However, the phase differences 
of the SRR oscillators in the array with respect to the 
driving field are generally different for each oscillator, as 
can be observed in Fig. 14, where the time-evolution of 
qi — qi is followed for approximately two periods Tb of 
the DDB oscillation. Note also that the time-evolution 
seems practically sinusoidal (harmonic), that may not be 
necessarily true for DDBs obtained with some other pa- 
rameter set. 

Dissipative Breathers hy frequency chirping.- For a fre- 
quency gapped linear spectrum, some of the linear modes 
become unstable at large amplitude. If the curvature of 
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FIG. 13: Dissipative breather profiles at maximum ampUtude 
for several values of the coupling parameter as shown on the 
figure which are localized at the surface (at n = 1), along with 
an almost uniform solution for A = 0.072 just above the value 
of \max for this particular parameter set. The parameters are 
n = 0.92, X = +1/6, 7 = 0.01, eo = 0.04, and 5 = 0.8. 




X 



FIG. 14: (Color online) Time-dependence of qi (black solid 
curve), 52 (red dotted curves), (green short-dashed curve), 
and qi (blue dashed curves), for a surface breather localized 
at n = 1 and A = 0.03. The other parameters are the same 
with those in the caption of Fig. 14. 



the dispersion curve in the region of that mode is neg- 
ative and the lattice potential is hard then, the large 
amplitude mode becomes unstable with respect to for- 
mation of a DB in the gap above the linear spectrum 
13; 20]. Below we exploit MI in order to generate spon- 
taneously DDBs in the binary array. The procedure that 
is followed is shortly described below [Ij, [20]. For the 
parameters in the captions of Fig. 15 and 16, the top of 
the upper linear band is located at ~ 1.42 where the 
curvature is negative. Moreover, the SRRs are subjected 
to on-site potentials that are hard (for % < 0) . The (large 
amplitude) driver is initiated with its frequency just be- 
low il and is then chirped with time to produce enough 
vibrational amplitude to induce MI of the uniform mode, 
which then leads to spontaneous DDB generation. At the 
end of the frequency chirping phase, the driver frequency 
is well above f2, and only supplies energy into the formed 
DDB(s). During that phase, a large number of DDBs 
may be generated, which can move and collide and even- 
tually coalesce into a small number of high amplitude 
DDBs that are frequency locked to the driver and, be- 
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parameter regime these breathers may be left-handed in 
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FIG. 15: Contours of the energy density 7i„ on the t — n 
plane for a strongly driven binary magnetic metamaterial, 
with 6 ^ 2, n = 1.42, x = -1/6, 7 = 0.001, A = 0.05, and 
eo = 3.0. 

cause of that, they are trapped at particular SRRs. After 
that, the driver frequency is kept constant and the high 
amplitude DDBs (and even some low-amplitude ones) 
continue to receive energy falling into a stationary state. 
When the driver is switched off all DDBs die out in a 
short time interval. 

In Figs. 15 and 16, the contours of the energy density 
Hn on the r — n plane identify the evolution of the DDBs 
formed by that procedure. The chirping phase lasts for 
2000 To time units (Tq — 2n/Q), where the frequency 
varies hnearly from Q, = 0.997f} to = imOfl. The 
driver is subsequently kept at constant frequency f2 / un- 
til it is switched off after another 2000 Tq time units. 
Figs. 15 and 16 correspond to the regions of the binary 
MMs where several DDBs have survived after the chirp- 
ing phase. In Fig. 15 we clearly observe a high amplitude 
DDB at rt = 251, along with some other DDBs of consid- 
erably lower amplitude, that survive until the end of the 
constant frequency phase. In Fig. 16, where the binary 
MM is driven not as strongly as that in Fig. 15, we ob- 
serve one relatively low-amplitude DDB which however 
survives until the end of the constant frequency phase. 
It is possible that the procedure described above, which 
relies on the MI of the large amplitude linear modes, can 
be used for the generation of DDBs in other magnetoin- 
ductive systems as well d^Ell) in their binary versions. 

VI. CONCLUSION 

We presented detailed analysis for induced nonlinear 
localization in binary nonlinear magnetic metamaterials. 
The systems we analyzed are one dimensional and con- 
sist of two types of SRRs; this configuration leads to 
a linearized magnetoinductive system with two optical 
bands separated by a gap. When nonlinearity is also 
taken into account nonlinear localized modes of discrete 
breather type may be generated in the gaps. As in the 
case where all units are identical [2^,[3^ depending on the 



FIG. 16: Contours of the energy density Tin on the t — n 
plane for a strongly driven binary magnetic metamaterial, 
with S = 2, n ^ 1.42, X = -1/6, 7 = 0.001, A = 0.05, and 
eo = 2.0. 



a right-handed background or the opposite. We focussed 
on both the Hamiltonian as well as the dissipative case; 
the latter is clearly the most interesting physically since it 
corresponds to true propagation of waves in the medium. 
We used two approaches, one based on the rotating wave 
approximation while the second on exact numerics us- 
ing the breather analysis from the anticontinuous limit. 
The comparison of the two showed that the RWA is in 
most cases a good approximation for a relatively accurate 
breather construction. 

In the Hamiltonian case we found two types of 
breathers with even or odd local symmetry depending 
on the sign of the product of the coupling parameter and 
the nonlinearity parameter. Both types may exist in the 
bulk but also in the boundary of the chain; the latter 
form surface breathers. A similar situation occurs also 
in the dissipative case where depending on the number 
of attractors of the single driven nonlinear oscillator sys- 
tem we have different type of dissipative breathers. Both 
bulk and surface breathers appear with corresponding 
symmetries. 

The binary structure of the lattice allows for genera- 
tion of breathers through direct external induction. This 
is accomplished through frequency chirping to the de- 
sired frequency. In the process of frequency modulation 
induction of plane wave instability occurs that leads to 
breather generation. These breathers move around in 
the lattice, collide, some decay and eventually a single 
"large" breather is left in the metamaterial. This method 
of breather generation is direct and may be used for ex- 
perimental breather investigation in metamaterials. 



Acknowledgments 

One of us (MIM) acknowledges support from Fondecyt 
Grant 1080374 of Chile. 



10 



[1] A. J. Sievers and S. Takeno, Phys. rev. Lett. 61, 970 
(1988). 

[2] For a recent review see S. Flach, A.V. Gorbach, Phys. 

Rep. 467, 1 (2008). 
[3] R. S. MacKay and S. Aubry, Nonlinearity 7, 1623 (1994). 
[4] S. Aubry, Physica D 103, 201 (1997). 
[5] J. L. Man'n and S. Aubry, Nonlinearity 9, 1501 (1996). 
[6] J. L. Man'n, F. Falo, P. J. Martinez, and L. M. Floria, 

Phys. Rev. E 63, 066603 (2001). 
[7] D. Zueco, P. J. Martinez, L. M. Floria, and F. Falo, Phys. 

Rev. E 71, 036613 (2005). 
[8] H. S. Eisenberg, Y. Silberberg, R. Morandotti, A. R. 

Boyd, and J. S. Aitchison, Phys. Rev. Lett. 81, 3383 

(1998). 

[9] B. I. Swanson, J. A. Brozik, S. P. Love, G. F. Strouse, A. 
P. Shreve, A. R. Bishop, W.-Z. Wang, and M. L Salkola, 
Phys. Rev. Lett. 82, 3288 (1999). 
[10] F. M. Russell, and J. C. Eilbeck, Europhys. Lett. 78, 
10004 (2007). 

[11] U. T. Schwarz, L. Q. English, and A. J. Sievers, Phys. 

Rev. Lett. 83, 223 (1999). 
[12] E. Trias, J. J. Mazo, and T. P. Orlando, Phys. Rev. Lett. 

84, 741 (2000). 

[13] M. Sato, B. E. Hubbard, A. J. Sievers, B. Ilic, D. A. 
Czaplewski, and H. G. Graighead, Phys. Rev. Lett. 90, 
044102 (2003). 

[14] M. Sato and A. J. Sievers, Phys. Rev. Lett. 98, 214101 
(2007). 

[15] G. P. Tsironis and S. Aubry, Phys. Rev. Lett. 77, 5225 
(1996). 

[16] K. O . Rasmussen, S. Aubry, A. R. Bishop and G. P. 

Tsironis, Eur. Jour. Phys. B 15, 169 (2000). 
[17] T. Dauxois and M. Peyrard, Phys. Rev. Lett. 70, 3935 

(1993). 

[18] D. Hennig, L. Schimansky-Geier, and P. Hanggi, Euro- 
phys. Lett. 78, 20002 (2007). 

[19] D. Hennig, S. Fugmann, L. Schimansky-Geier, and P. 
Hanggi, Acta Physica Polonica B 39, 1125 (2008). 

[20] M. Sato, B.E. Hubbard, A.J. Sievers, Rev. Mod. Phys. 
78 (2006) 137. 

[21] M. Gorkunov, M. Lapirie, E. Shamonina, and K. H. 

Ringhofer, Eur. Phys. J. B 28, 263 (2002). 
[22] S. Linden, C. Enkrich, G. Dolling, M. W. Klein, J. Zhou, 

T. Koschny, C. M. Soukoulis, S. Burger, F. Schmidt, and 

M. Wegener, IEEE J. Selec. Top. Quant. Electron. 12, 

1097 (2006). 

[23] T. J. Yen, W. J. Padilla, N. Fang, D. C. Vier, D. R. 



Smith, J. B. Pendry, D. N. Basov, and X. Zhang, Science 
303, 1494 (2004). 
[24] N. Katsarakis, G. Konstantinidis, A. Kostopoulos, R. S. 
Penciu, T. F. Gundogdu, M. Kafesaki, E. N. Economou, 
Th. Koschny, and C. M. Soukoulis, Opt. Lett. 30, 1348 

(2005) . 

[25] T. H. Hand and S. A. Cummer, J. Appl. Phys. 103, 
066105 (2008). 

[26] D. A. Powell, I. V. Shadrivov, Yu. S. Kivshar, and M. V. 

Gorkunov, Appl. Phys. Lett. 91, 144107 (2007). 
[27] I. V. Shadrivov, S. K. Morrison, and Yu. S. Kivshar, Opt. 

Express 14, 9344 (2006). 
[28] I.V. Shadrivov, A.B. Kozyrev, D. van der Weide, Yu.S. 

Kivshar, Appl. Phys. Lett. 93, 161903 (2008). 
[29] N. Lazarides, M. Eleftheriou, and G. P. Tsironis, Phys. 

Rev. Lett. 97, 157406 (2006). 
[30] M. Eleftheriou, N. Lazarides, and G. P. Tsironis, Phys. 

Rev. E 77, 036608 (2008). 
[31] N. Lazarides, G. P. Tsironis, and Yu. S. Kivshar, Phys. 

Rev. E 77, 036608 (2008). 
[32] M. Eleftheriou, N. Lazarides, G. P. Tsironis, and Yu. S. 

Kivshar, arXiv:0903.2130vl [cond-mat.mtrl-sci] (subm. 

to Phys. Rev. E). 
[33] I. V. Shadrivov, A. A. Zharov, N. A. Zharova, and Yu. 

S. Kivshar, Photonics Nanostruct. Fundam. Appl. 4, 69 

(2006) . 

[34] I. Kourakis, N. Lazarides, and G. P. Tsironis, Phys. Rev. 

E 75, 067601 (2007). 
[35] M. L Molina, R. A. Vicencio, and Yu. S. Kivshar, Opt. 

Lett. 31, 1693 (2006). 
[36] Yu. S. Kivshar and M. L Molina, Wave Motion 45, 59 

(2007) . 

[37] M. V. Gorkunov, L V. Shadrivov, and Yu. S. Kivshar, 

Appl. Phys. Lett. 88, 071912 (2006). 
[38] E. Shamonina, V. A. Kalinin, K. H. Ringhofer, and L. 

Solymar, J. Appl. Phys. 92, 6252 (2002). 
[39] I. V. Shadrivov, A. N. Reznik, and Yu. S. Kivshar, Phys- 
ica B 394 180 (2007). 
[40] A. A. Zharov, I. V. Shadrivov, and Y. S. Kivshar, Phys. 

Rev. Lett. 91, 037401 (2003). 
[41] M. I. Molina, L L. Garanovich, A. A. Sukhorukov, and 

Yu. S. Kivshar, Opt. Lett. 31, 2332 (2006). 
[42] N. Lazarides, G. P. Tsironis, and M. Eleftheriou, Nonlin. 

Phen. Compl. Syst. 11, 250 (2008). 
[43] G. P. Tsironis, N. Lazarides, and M. Eleftheriou, PIERS 

Online 5, 26 (2009). 



